# Clear memory
rm(list = ls())

# Load packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  sf, maps, data.table, janitor, magrittr, fixest,
  broom, tidyverse, tidylog
)
options("tidylog.display" = NULL)
`%notin%` <- Negate(`%in%`)
transf <- function(x) (exp(x - 1) - 1) * 100

# # ===========================
# ##### value of adaptation
# # ===========================

welfare_base <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:unemp_amenity, ~ transf(.x))) %>%
  pull(total, consumption_total)

welfare_all <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-off_mobility-off_mktacc-off_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:unemp_amenity, ~ transf(.x))) %>%
  pull(total)

welfare_labor <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-off_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:unemp_amenity, ~ transf(.x))) %>%
  pull(total)

welfare_trade <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-off_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:unemp_amenity, ~ transf(.x))) %>%
  pull(total, consumption_total)

welfare_phys <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-off_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(across(total:unemp_amenity, ~ transf(.x))) %>%
  pull(total)

# # ===========================
# ##### 25 counties account for welfare
# # ===========================

workers <- fread("data/payroll_cf.csv") %>% as_tibble()

# total number of workers in billions
workers_total <- sum(workers$employment, na.rm = T) / 1e9

welfare_total <- fread("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  select(fips, industry, welfare, cf_labor, cf_wage, real_wage, cf_c_price, amenities) %>%
  # consumption equivalent change in welfare (100%) * base consumption * (1%/100%) * (10000 $ / 1$)
  mutate(total_wage_gain = transf(welfare) * cf_wage / cf_c_price / 100 * 1e4) %>%
  group_by(fips) |>
  summarise(total_income_gain = sum(cf_labor * total_wage_gain * workers_total, na.rm = TRUE))

# 25 counties account for half of the $41 billion in welfare

#View(welfare_total |> arrange(desc(total_income_gain)) |> mutate(cumtotal = cumsum(total_income_gain)))
sum(welfare_total$total_income_gain)
print(n = 25, head(welfare_total |> arrange(desc(total_income_gain)) |> mutate(cumtotal = cumsum(total_income_gain)), n = 25))
# 25 account for half of $41 billion 

# ===========================
##### Where manufacturing moves
# ===========================

# change relative to no nonattainment
employment_change <- fread("data/counterfactual/welfare_total_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() |>
  select(m_emp_hat:u_emp_percent)
employment_change
employment_change$m_emp_hat

percent_to_unemp <- employment_change$u_emp_hat / employment_change$o_emp_hat

percent_to_unemp

# ===========================
##### Where manufacturing moves under first-best
# ===========================

# change in first-best relative to 1997 nonattainment
employment_change <- fread("data/counterfactual/welfare_total_opt-opt_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() |>
  select(m_emp_hat:u_emp_percent)
employment_change
employment_change$m_emp_hat

percent_to_unemp <- employment_change$u_emp_hat / employment_change$o_emp_hat

percent_to_unemp


# ===========================
##### Percent of population
##### better off with first-best
# ===========================

first_best_pop_improvement <- read.csv("data/counterfactual/welfare_opt-opt_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  as_tibble() %>%
  # opt models are solved in the opposite relative way so we want welfare > 1 for them to
  # be worse off under the 1997 NAAQS
  mutate(worse_off = welfare > 1) %>%
  group_by(worse_off) %>%
  summarise(base_labor = sum(base_labor))
first_best_pop_improvement

# ===========================
##### Emissions Leakage
# ===========================

fips <- fread("data/simulation_fips_crosswalk.csv")

non_attainment <- fread("data/productivity_panel.csv") %>%
  filter(year >= 1986 & year <= 1997) %>%
  group_by(fips, industry) %>%
  mutate(total_nonattainment = sum(nonattainment, na.rm = T)) %>%
  ungroup() %>%
  mutate(nonattainment_1991 = nonattainment > 0 & year == 1991) %>%
  mutate(nonattainment_1997 = nonattainment > 0 & year == 1997) %>%
  group_by(fips) %>%
  mutate(nonattainment_1991 = max(nonattainment_1991)) %>%
  mutate(ever_nonattainment = max(nonattainment_1997) > 0) %>%
  ungroup() %>%
  filter(year == 1997) %>%
  mutate(fips = str_pad(fips, 5, "left", "0")) %>%
  distinct(fips, nonattainment_1997)

county_fips <- county.fips
counties_sdf <- st_as_sf(maps::map("county", plot = FALSE, fill = TRUE)) %>%
  rename(polyname = ID) %>%
  full_join(county_fips) %>%
  mutate(fips = ifelse(fips == 46113, 46102, fips)) %>%
  mutate(fips = ifelse(polyname == "florida,okaloosa", 12091, fips)) %>%
  mutate(fips = ifelse(polyname == "virginia,accomack", 51001, fips)) %>%
  mutate(fips = ifelse(polyname == "texas,galveston", 48167, fips)) %>%
  mutate(fips = ifelse(polyname == "louisiana,st martin", 22099, fips)) %>%
  mutate(fips = ifelse(polyname == "north carolina,currituck", 37053, fips)) %>%
  mutate(fips = ifelse(polyname == "washington,pierce", 53053, fips)) %>%
  mutate(fips = ifelse(polyname == "washington,san juan", 53055, fips)) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  st_transform(5070)

welfare_base <- fread("data/counterfactual/welfare_opt-base_prod-0.0_amen-on_abatement-on_transport-on_mobility-on_mktacc-on_cong-off_heterog-off_info-on_theta-4.0_alpha-0.274_gamma-0.481_iota-1.0_xi_mult-1.0_vsl-cons.csv") %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  left_join(counties_sdf) %>%
  left_join(non_attainment) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  mutate(
    emissions_nh3_change = emissions_nh3 * cf_emissions_nh3 - cf_emissions_nh3,
    emissions_nox_change = emissions_nox * cf_emissions_nox - cf_emissions_nox,
    emissions_pm25_change = emissions_pm25 * cf_emissions_pm25 - cf_emissions_pm25,
    emissions_so2_change = emissions_so2 * cf_emissions_so2 - cf_emissions_so2,
    emissions_voc_change = emissions_voc * cf_emissions_voc - cf_emissions_voc
  )

leakage <- welfare_base %>%
  group_by(nonattainment_1997) %>%
  summarise(across(emissions_nh3_change:emissions_voc_change, ~ sum(.x))) %>%
  mutate(across(emissions_nh3_change:emissions_voc_change, ~ .x / sum(.x), .names = "{.col}_share_of_decline"))
leakage
fips <- fread("raw/modified-ap3/fips.csv")$V1

nh3_damages <- fread("data/modified-ap3/NH3_md_L_1997_prime_aged.csv") %>%
  as_tibble()
nox_damages <- fread("data/modified-ap3/NOX_md_L_1997_prime_aged.csv") %>%
  as_tibble()
pm25_damages <- fread("data/modified-ap3/PM25_md_L_1997_prime_aged.csv") %>%
  as_tibble()
so2_damages <- fread("data/modified-ap3/SO2_md_L_1997_prime_aged.csv") %>%
  as_tibble()
voc_damages <- fread("data/modified-ap3/VOC_md_L_1997_prime_aged.csv") %>%
  as_tibble()

names(nh3_damages) <- as.character(fips)
names(nox_damages) <- as.character(fips)
names(pm25_damages) <- as.character(fips)
names(so2_damages) <- as.character(fips)
names(voc_damages) <- as.character(fips)

nh3_damages <- nh3_damages %>%
  mutate(fips = fips) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  relocate(fips) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  mutate(total_damage = rowSums(across(where(is.numeric)))) %>%
  select(fips, nh3_damage = total_damage)

nox_damages <- nox_damages %>%
  mutate(fips = fips) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  relocate(fips) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  mutate(total_damage = rowSums(across(where(is.numeric)))) %>%
  select(fips, nox_damage = total_damage)

pm25_damages <- pm25_damages %>%
  mutate(fips = fips) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  relocate(fips) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  mutate(total_damage = rowSums(across(where(is.numeric)))) %>%
  select(fips, pm25_damage = total_damage)

so2_damages <- so2_damages %>%
  mutate(fips = fips) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  relocate(fips) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  mutate(total_damage = rowSums(across(where(is.numeric)))) %>%
  select(fips, so2_damage = total_damage)

voc_damages <- voc_damages %>%
  mutate(fips = fips) %>%
  mutate(fips = str_pad(fips, 5, "left", pad = "0")) %>%
  relocate(fips) %>%
  mutate(fips = ifelse(fips == "12025", "12086", fips)) %>%
  mutate(total_damage = rowSums(across(where(is.numeric)))) %>%
  select(fips, voc_damage = total_damage)

damage_leakage <- welfare_base %>%
  select(fips, nonattainment_1997, emissions_nh3_change:emissions_voc_change) %>%
  left_join(nh3_damages) %>%
  left_join(nox_damages) %>%
  left_join(pm25_damages) %>%
  left_join(so2_damages) %>%
  left_join(voc_damages) %>%
  mutate(
    nh3_damage = nh3_damage * emissions_nh3_change,
    nox_damage = nox_damage * emissions_nox_change,
    pm25_damage = pm25_damage * emissions_pm25_change,
    so2_damage = so2_damage * emissions_so2_change,
    voc_damage = voc_damage * emissions_voc_change
  ) %>%
  replace_na(list(nonattainment_1997 = FALSE)) %>%
  group_by(nonattainment_1997) %>%
  summarise(across(nh3_damage:voc_damage, ~ sum(.x)))

damage_leakage
